// -*- coding: utf-8 -*-

#include <iostream>
#include <string>

#include "Modeler.h"

using namespace std;
using namespace elisa;

/**
 * Lorenz方程式と球
 */

class Lorenz : public Modeler
{
private:
	static Interval Time;
	static VNInterval X1, X2, X3;

	ContVar x1, x2, x3;

public:
	Lorenz()
		: Modeler(Time),
		  x1(X1, "x1"),
		  x2(X2, "x2"),
		  x3(X3, "x3")
	{
		//elisa::RealVar::SetPrecision(1.0e-11);
		elisa::RealVar::SetPrecision(1.0e-1);
		t.setPrecision(1.0e-2);
		t.setDomain(new Interval(Time.inf(), 20.0));

		addVar(x1);
		addVar(x2);
		addVar(x3);
	}

	// 連続制約の定義
	template<typename T>
	static void flow(int n, T *dv, const T *v, T t, void *param)
	{
		const T Sigma(10);
		const T Rho  (28);
		const T Beta (T(8)/T(3));

		dv[0] = Sigma * (v[1]-v[0]);
		dv[1] = v[0] * (Rho-v[2]) - v[1];
		dv[2] = v[0]*v[1] - Beta*v[2];
	}

	// ガード制約の定義
	template<typename T> 
	static T guard(const std::vector<T>& v, const T& t) {
		return (pow(v[0],2) + pow(v[1],2) + pow((v[2]-T(28)),2)) - 700;
	}

	bool reset(bool first)
	{
		return first ? true : false;
	}
	
	void print(std::ostream& out) const
	{
		cout << '\t' << t.getName()  << " = " << t.domain()  << endl;
		cout << '\t' << x1.getName() << " = " << x1.domain() << endl;
		cout << '\t' << x2.getName() << " = " << x2.domain() << endl;
	}
};

MODEL(Lorenz)


Interval Lorenz::Time(Interval::Zero);
VNInterval Lorenz::X1(15);
VNInterval Lorenz::X2(15);
VNInterval Lorenz::X3(36);
